A network-based modeling framework reveals the core signal transduction network underlying high carbon dioxide-induced stomatal closure in guard cells

Stomata are pores on plant aerial surfaces, each bordered by a pair of guard cells. They control gas exchange vital for plant survival. Understanding how guard cells respond to environmental signals such as atmospheric carbon dioxide (CO2) levels is not only insightful to fundamental biology but also relevant to real-world issues of crop productivity under global climate change. In the past decade, multiple important signaling elements for stomatal closure induced by elevated CO2 have been identified. Yet, there is no comprehensive understanding of high CO2-induced stomatal closure. In this work, we assemble a cellular signaling network underlying high CO2-induced stomatal closure by integrating evidence from a comprehensive literature analysis. We further construct a Boolean dynamic model of the network, which allows in silico simulation of the stomatal closure response to high CO2 in wild-type Arabidopsis thaliana plants and in cases of pharmacological or genetic manipulation of network nodes. Our model has a 91% accuracy in capturing known experimental observations. We perform network-based logical analysis and reveal a feedback core of the network, which dictates cellular decisions in closure response to high CO2. Based on these analyses, we predict and experimentally confirm that applying nitric oxide (NO) induces stomatal closure in ambient CO2 and causes hypersensitivity to elevated CO2. Moreover, we predict a negative regulatory relationship between NO and the protein phosphatase ABI2 and find experimentally that NO inhibits ABI2 phosphatase activity. The experimental validation of these model predictions demonstrates the effectiveness of network-based modeling and highlights the decision-making role of the feedback core of the network in signal transduction. We further explore the model’s potential in predicting targets of signaling elements not yet connected to the CO2 network. Our combination of network science, in silico model simulation, and experimental assays demonstrates an effective interdisciplinary approach to understanding system-level biology.


Introduction
Stomata are microscopic pores on plant leaf surfaces that control gas exchange such as the uptake of CO 2 and the release of water vapor and photosynthetically produced O 2 .Each stomate is bordered by 2 guard cells, whose shape change in response to the external and internal environment controls the stomate's opening or closure.Because plants inevitably lose water via transpiration through the stomata, the size of the stomatal pores is tightly regulated to balance the competing needs of CO 2 supply for photosynthesis and plant hydration.Stomata open in response to environmental signals such as light and low atmospheric CO 2 levels and close in response to darkness, high CO 2 , drought, and low humidity [1].As elucidated in detail for the stress hormone abscisic acid (ABA) [2][3][4], stimulus perception activates an intracellular signal transduction network, which leads to Ca 2+ influx and anion efflux, with consequent membrane depolarization promoting K + efflux.The outflow of ions increases cellular water potential, creating a driving force for water efflux through aquaporins, which finally causes guard cell deflation and stomatal closure.The extent to which the ABA signal transduction network is shared by other closure signals such as darkness and high CO 2 has been studied and debated [2,5].
Understanding how stomata respond to high CO 2 is a fundamental question for plant biology.Since stomata control gas exchange vital for plant water retention and photosynthesis, stomatal response to abiotic stress is an essential process for plant survival.Furthermore, understanding plant response to high CO 2 also has imperative practical significance: over the past 60 years, CO 2 concentration in the atmosphere has increased by 32% (NOAA at Mauna Loa [6]).The greenhouse effect induced by increased CO 2 has increased global temperature, which exacerbates drought, compounding the stress.Plant stomatal responses to CO 2 are therefore highly relevant to real-world issues of crop productivity under the adversity posed by climate change.

Current knowledge of high CO 2 -induced stomatal closure in guard cells
Stomatal movement in response to high CO 2 concentration is a complex process in which numerous cellular elements have been implicated by prior studies.These elements interact with and regulate each other via chemical reactions, protein interactions, and physical processes, forming a complex cell signaling network that propagates the signal from CO 2 to stomatal closure.Atmospheric CO 2 enters guard cells through CO 2 porins [7].The CO 2 -binding carbonic anhydrase enzymes CA1, CA4, and CA7 [8][9][10] catalyze the reversible chemical reaction CO 2 + H 2 O $ HCO 3 -+ H + , and these enzymes are a demonstrated component of the CO 2 signal transduction mechanism.HT1 (HIGH LEAF TEMPERATURE 1) is an important CO 2 signaling element that inhibits the activation of the main anion channel SLAC1 (SLOW ANION CHANNEL-ASSOCIATED 1).Multiple mutants of HT1 have been shown to result in insensitivity to high CO 2 [11][12][13].Mitogen-activated protein kinases 4 and 12 (MPK4/ MPK12) bind to HT1 and inhibit its kinase activity [13,14].Bicarbonate enhances and stabilizes the binding of MPK4/MPK12 to HT1 [15,16], resulting in loss of HT1 activation of CBC1 kinase, which itself is an inhibitor of closure [15].Arabidopsis RESISTANT TO HIGH CO 2 (RHC1), a MATE-type transporter, may also link elevated CO 2 concentration to repression of HT1, although this has been debated [17,18].Disruption of the second messengers Ca 2+ cyt , ROS (reactive oxygen species), or NO (nitric oxide) impairs high CO 2 -induced closure [19][20][21][22][23][24], implicating them in CO 2 signaling.Despite the recent identification of key elements and pathways in the high CO 2 response, further elucidation is needed of how the known elements and pathways interact with each other to induce stomatal closure.A systems biology approach integrating existing evidence into a signaling network can identify gaps in knowledge and suggest potential ways to elucidate all the mechanisms that contribute to stomatal closure.

Prior modeling of guard cell signaling networks
In the past, our team has developed multiple networks and network-based models to characterize the guard cell signaling system.In 2006, Li and colleagues constructed a Boolean model of ABA-induced stomatal closure [25].In 2014, Sun and colleagues developed a multilevel discrete dynamic model focusing on light-induced stomatal opening [26].In 2017, Albert and colleagues made major updates to the ABA-induced stomatal closure model, including many more signaling elements discovered in recent years [3].These network-based models revealed mechanistic insights and made novel predictions that were validated by follow-up experiments.Furthermore, we have developed novel theoretical concepts such as the concept of the stable motif, which forms the foundation of computational methods to determine the repertoire of long-time behaviors of biological systems and to identify nodes that can drive the system into a desired state [27][28][29][30].
Karanam and colleagues recently developed a graphical user interface to easily analyze Boolean models [31] and applied it to integrate a CO 2 signaling pathway into our previous ABA-induced closure model [3], as an approach to describe CO 2 signaling.While it is a good demonstration of their methodology, this preliminary model is not a rigorous representation of CO 2 signaling in that: (1) it lacks a comprehensive literature review to identify elements that are differentially regulated in high CO 2 signaling as compared to ABA signaling; and (2) it assumes that high CO 2 is sufficient to directly activate water efflux through the membrane, as well as Ca 2+ inflow, which together are sufficient to drive stomatal closure: these assumptions short-circuit all other CO 2 -related nodes, i.e., render the rest of their network superfluous.Because of the assumed direct connection to water efflux, this model would fail to recapitulate, for example, the experimental observations of impaired high CO 2 -induced closure for ht1 [12] and mpk4/mpk12 [14] mutants.A more rigorous network construction approach is required to produce a more accurate representation and simulation of high CO 2 -induced closure.
In this work, we compile and integrate an extensive body of literature evidence to assemble a guard cell signal transduction network that characterizes how high CO 2 induces stomatal closure through the decrease in turgor and volume of guard cell pairs.To further capture the dynamics involved in the stomatal closure process, we develop a dynamic model of the signaling network.We rigorously evaluate the accuracy of the model against known experimental observations, obtaining a high (91%) accuracy.We perform analyses including network connectivity, model simulation, and stable motif analysis, which reveal a feedback core that dictates the cellular decisions in the CO 2 response, and the roles of key nodes in the process.We then demonstrate the predictive power of our network model by predicting the outcomes of numerous perturbations (e.g., node knockouts) that can potentially influence stomatal response; such predictions are valuable for the prioritization of wet bench experiments.Finally, we experimentally test and validate 3 predictions related to the role of NO in stomatal closure.A flowchart of our work is shown below (Fig 1).Our modeling framework captures systemlevel biological processes and provides insights into their logic and key underlying mechanisms.

Construction and analysis of the signaling network of high CO 2 -induced stomatal closure
We conducted a comprehensive literature review, covering more than 160 publications relating to stomatal closure or to elements of guard cell CO 2 signaling, and manually curated the interaction information from these publications (S1 Table ), aiming to identify and integrate the mechanisms that lead to high CO 2 -induced rapid stomatal closure.For this purpose, guard cell development and long-term responses based on gene regulation are out of the scope of our model.Further, not all guard cell signaling elements are equally relevant to high CO 2 -induced closure.For example, cytosolic pH increases in response to ABA, but does not increase under high CO 2 [9,32,33].Another example is the kinase OST1 (Open Stomata1/SnRK2.6),known to phosphorylate SLAC1 as a key contributor to the ion flow that leads to stomatal closure induced by ABA.Although ost1 mutants show strongly impaired stomatal closure under both ABA and high CO 2 , recent studies showed that unlike ABA, high CO 2 does not increase OST1 kinase activity in guard cells [34,35].
To focus on CO 2 signaling, we divide the evidence into 5 categories based on relevance: (1) elements and their interactions that are validated to be involved in CO 2 signaling; (2) generic interactions that can be assumed to be involved in CO 2 signaling, e.g., chemical reactions whose reactants are reasonably expected to be present in the guard cell; (3) evidence that shows an element is not involved in CO 2 signaling; (4) evidence that supports an element's role in stomatal closure in response to another signal, but the element's involvement in high CO 2 -induced closure has not been investigated; and (5) evidence that supports an element's Flow chart of the process of network construction and analysis.We construct a signaling network based on literature review, resulting in a network of 47 nodes and 95 edges.We then construct a Boolean dynamic model, which adds node states and regulatory logic to the network.We assess the accuracy of our network model by comparing its results with known experimental observations, finding that the model simulations successfully recapitulate 91% of experimental observations.We also perform multiple network-based analyses on the model.Finally, we experimentally validate selected model predictions.
https://doi.org/10.1371/journal.pbio.3002592.g001role in high CO 2 -induced closure but this element has no or insufficient connections to elements in categories (1) and ( 2).The resulting list of observations and their categorization is provided in S1 Table .We use category ( 1) and ( 2) evidence to construct the high CO 2 signaling network and exclude category (3) and ( 4) evidence.Category (5) evidence is also necessarily excluded, as it is not possible to connect such nodes to the network.Such categorization is also helpful to wet bench investigators; e.g., category 4 and 5 elements would be of particular interest to investigate further in the future.
To represent signaling elements with different functions, we deploy a multi-node representation for the OST1, MPK, and SLAC1 proteins, separating the different (e.g., CO 2 -dependent and independent) mechanisms for their activation (see Fig 2C for an example, and see MPKs are known to inhibit HT1 by a physical interaction induced by bicarbonate and independent of the MPKs' kinase activity [15,16], and also are necessary in Ca 2+ -induced SLAC1 activation [37].https://doi.org/10.1371/journal.pbio.3002592.g002

S2 Table and S1 Text for details).
As the follow-up analyses require that each edge of the network points to a node, during network construction we redirected each evidence of influence on an edge to an equivalent influence on the target of the edge [36].
Biological description of the high CO 2 -induced stomatal closure network and its overlap with the ABA-induced closure network.The signaling network corresponding to high CO 2 -induced closure contains 47 nodes and 95 edges (Fig 2A ).A summary table of all edges is presented in Table 1.The full list of nodes and edges, including the full node names, is provided in S2 Table .We mark nodes of the network (Fig 2A) in 3 colors, to reflect their specificity to high CO 2 -induced closure.Yellow nodes represent validated elements in CO 2 signaling, i.e., category (1); white nodes represent assumed elements based on generic evidence, i.e., category (2).For example, the white nodes 8-nitro-cGMP, ADPRc, and cADPR represent metabolites known to be produced in response to NO (a validated CO 2 element) [38] and are necessary to connect NO to the rest of the network.Although by default nodes in categories (3) to (5) are not included in the network, "OST1 activated" is a member of a node pair (it represents increased kinase activity of OST1), thus we include it and mark it in purple to indicate that it belongs to category (4): it is not involved in high CO 2 -induced closure but is important under ABA-induced closure.
The edges in the network are also marked to distinguish their properties: solid edges indicate a direct interaction/regulation, dashed edges indicate indirect regulation, and dotted edges indicate assumed regulation.An example of indirect regulation is the edge from ROS to the H + ATPase, which is dashed.An example of an edge referring to assumed regulation is the edge from HCO 3 -to SLAC1 CO 2 , which is supported by a specific SLAC1 residue that is required for full CO 2 signaling and is structurally predicted to interact with HCO 3 - [60], but this predicted physical interaction has not been experimentally verified.Like this example, all assumed edges in the model are well-supported.Further justification of the assumed edges can be found in S1 Text.
The high CO 2 -induced stomatal closure network is known to share nodes and pathways with the ABA-induced closure network, but it also has its unique signaling elements.Furthermore, through the literature review process, we also found multiple elements that participate in ABA-induced closure but are not yet evaluated in the context of high CO 2 -induced closure (also see S1 Table and S2 Text).We list in Table 2 the different categories for the combined set of the signaling elements (i.e., nodes) of this CO 2 network and of our 2017 ABA network [3].
The CO 2 signaling network consists of 3 distinct modules.Of the 47 nodes in the network, 8 (including high CO 2 ) are source nodes that are not regulated by any other node, and 1 node (Closure) is a sink node (output).The source nodes other than the signal "high CO 2 " mainly represent molecules needed as substrates for reactions incorporated in the network: CO 2 porins, CAs, GTP, MPKs, NAD + , NADPH, Nitrite.Structural and connectivity analysis shows that the CO 2 network is organized into 3 graph modules (called "connected components" in graph theoretical language), as shown in Fig 2B .First, 7 early CO 2 signaling nodes form an in-component, i.e., these nodes are not regulated by any nodes from other components in the network.This module mainly represents how CO 2 enters the guard cell and is converted to HCO 3 -.Second, a large strongly connected component (SCC) of 20 nodes forms the core of the signal transduction network."Strongly connected" means that for any node pair A, B, there is at least 1 regulatory path from A to B and at least 1 path from B to A. The feedback loops underlying the strong connectivity allows complex dynamical behavior that ultimately dictates the homeostasis of the system.The 2 edges between the in-component and the SCC represent CO 2 sensing mechanisms, e.g., the inhibitory effect of HCO 3 -on HT1 via the formation of the MPKs-HT1 complex or via RHC1.Third, 13 nodes, mainly representing ion transport mechanisms and the ion and water fluxes they mediate, as well as the output described in [60] forms a unique edge that connects the CO 2 in-component directly to the out-component.
Construction of the dynamic model.A network of signaling elements reveals the interactions and regulatory relationships between them, yet this alone is not sufficient to capture the dynamic aspects of signal transduction.In the stomatal closure process, the biological system represented by the network evolves from a state corresponding to open stomata (in the absence of a closure signal) to a state corresponding to closed stomata after sufficient time has elapsed in the presence of a closure signal.To capture the dynamic state transition of the CO 2 network, we perform dynamic modeling, where we describe each node of the network with a state variable that represents its abundance or activity.In this way, the propagation of the high CO 2 signal is captured in the model by a series of node state changes propagating through the network, eventually leading to stomatal closure.We use a Boolean model, i.e., we describe each node with binary states: ON or 1 for higher-than-threshold abundance or activity; OFF or 0 for lower-than-threshold abundance and activity.Despite its simplicity, Boolean dynamics has proved capable of capturing complex biological behaviors [74], e.g., cell differentiation as a consequence of multi-stability [75,76], or sustained oscillations such as the cell cycle [77].For each node, its state change is determined by the regulatory influences incident on it, which are captured by its Boolean regulatory function.The regulatory function uses the states of the node's regulators in the network as inputs and combines these inputs using logical operators such that the function correctly reflects experimental observations.For example, the node "HCO 3 -" represents bicarbonate concentration in the guard cell.Its input nodes are "CO 2 entry of cell" and "CAs" (the CO 2 -binding carbonic anhydrase enzymes CA1, CA4, and CA7).Since bicarbonate is produced from the reaction CO 2 + H 2 O !HCO 3 -+ H + , catalyzed by CAs, we know that both CO 2 entry into the cell and CAs are required for the increase of intracellular bicarbonate concentration, thus the inputs are connected with a logical "AND" operator, and the regulatory function of node "HCO 3 -" is set as "f HCO 3 À = CAs AND CO 2 entry of cell."The comprehensive list of regulatory functions and their justifications are provided in S1 Text.Each node's regulatory function determines the future state of the node as determined by the current state of its regulators.
With the Boolean regulatory functions defined for each node, we can perform simulations of the signal transduction process over the entire CO 2 network.We initialize the in silico stomata in a pre-stimulus, i.e., open state, with an ON state for closure-inhibiting nodes (namely ABI1, ABI2, CBC1/2, HAB1, HT1, and H + ATPase) and other nodes OFF (see S1 Text for the complete initial condition).We deploy a discrete-time simulation, with stochastic random order asynchronous update, and simulate the network model for 30 timesteps (see Methods).Each simulated trajectory of the model will eventually converge into a stabilized behavior of the system as a whole, called an attractor.The most intuitive attractor is a steady state, where all nodes stabilize in a fixed state.Alternatively, a system can also have an oscillatory attractor, where a fraction of nodes indefinitely oscillates between 0 (OFF) and 1 (ON).A system can have multiple attractors (called multi-stability), which may be reachable from the same initial state if the system's timing is stochastic.We obtain the consensus behavior of the system by performing 1,000 replicate simulations, i.e., 1,000 "in silico stomata" simulations starting in the same pre-stimulus state (see Methods), then calculating the percentage of node activation at each timestep as our metric, meaning the percentage of simulations in which the node (e.g., Closure) is in state 1 (ON).In this way, oscillations and multi-stability with different node states will display a percentage of node activation with a value between 0% and 100%.We demonstrate simulated high CO 2 -induced stomatal closure and lack of closure under ambient CO 2 in Fig 3B and show a summary of node states in the corresponding attractors in Fig 3C .In the simulated wild-type system (Fig 3B), the percentage of closure in response to high CO 2 increases over the timesteps from 0% to 100%.Despite the deliberate choice of stochasticity in the timing, all the simulations converge into the same attractor, in which the node Closure is ON.The row marked "high CO 2 " in Fig 3C shows the states of all nodes in this attractor, recapitulating the biological knowledge regarding the activity of nodes in stomata that closed in response to high CO 2 .We marked with "e" the node states that have been studied experimentally (see the entries in S1 Table in which Entity/Process A is High CO 2 and which report measurements of Entity/Process B).For example, the node "OST1 activated" is OFF in this closure attractor, meaning that OST1 does not become activated under high CO 2 , consistent with recent publications [34,35].The node Ca 2+ cyt oscillates between ON and OFF, capturing transient cytosolic Ca 2+ concentration increases that have been observed during CO 2 response [19,20].Node states for which no prior experiments exist represent novel predictions.For example, the model indicates that high CO 2 activates not only the CO 2 -specific "SLAC1 CO 2 " node but also the node "SLAC1 generic," as all of these node's regulators, including CPKs, are indirectly activated by high CO 2 .In the simulation curve corresponding to ambient CO 2 (represented as the inactivity of the node "high CO 2 "), the percentage of closure remains 0% for the entire duration of the simulation ( Validation of the network model.We assembled all experimental evidence related to high CO 2 -induced stomatal closure that was not used in model construction (see the entries with "V" for validation, in S1 Table ), and checked if the model simulations of the relevant settings yield results consistent with the experimental evidence.An important concept that can be insufficiently appreciated by experimentalists is that all the pathway-level effects (e.g., the effect of a node manipulation on high CO 2 -induced stomatal closure) are emergent properties of the system that arise from multiple interactions; thus, the fact that the regulatory functions of each node are informed by experimental evidence is not sufficient to guarantee consistency with pathway-level observations.The summary of validation and consistency is presented in Table 3.In a total of 47 settings, there are 43 cases of consistency between the experimental and simulation results and 4 inconsistent cases, yielding an accuracy of 43/47 = 91%.The high validation accuracy shows that our model correctly captures most known signaling behaviors in CO 2 -induced stomatal closure.Note that a few comparisons could be consistent by construction if the observed regulator node directly regulates the target node in the model.When assessing the validity of our model, we exclude this type of trivial consistency.A more detailed validation table including the trivially consistent cases, a description of each experimental observation, and explanations are provided in S3 Table.

Network-based analysis identifies the dynamic repertoire and decisionmaking mechanisms of the CO 2 signaling system
In the previous section, we showed that our model can accurately reproduce existing experimental observations.Next, we explore the explanatory power of the model and its predictions about the biological system's response in conditions for which experimental observations are not yet available.We start by determining the full attractor repertoire of the model.The attractor repertoire of the model is important because attractors correspond to biological phenotypes, e.g., open or closed stomata are represented by 2 different attractors in the CO 2 network model.Finding the attractor repertoire is the first step toward understanding how a biological system makes decisions on how to switch between different phenotypes.Understanding the decision-making behind attractor transitions of a biological system can reveal the key signal transduction mechanisms and suggest potential interventions to guide the system into a more desirable state.Finding attractors by simulations is difficult because the state space of a network model increases exponentially with its size, e.g., the CO 2 network model with 47 nodes has a total of 2 47 states.An alternative and efficient way to find the attractors of a system is to utilize the network topology via stable motif analysis, developed by the Albert group [86].Stable motifs are special feedback loops in a dynamic model that can self-sustain once their constituent nodes attain an associated state (see Methods).The important feature of stable motifs is that they are irreversible: after locking in they become independent of the rest of the system.This feature implies that the lock-in of a stable motif is like a decision of the system, and a series of stable motif lock-ins will lead to an attractor [28,30,86].Therefore, finding all stable motifs and their lock-in sequence allows comprehensive identification of attractors.Stable motif-based attractor identification was shown to be more efficient [29] than methods based on mapping state transitions [87][88][89], algebraic methods [90], and methods based on network decomposition [91].In addition, stable motif-based control outperforms other approaches in identifying interventions that can drive the system into a beneficial attractor [29].
Stable motifs identify the dynamic repertoire and decision-making of the system.Whether in ambient or high CO 2 , the unperturbed (wild type) system has a stable motif (see Fig 4A) that includes the activity of RbohD/F and ROS, the minimum functional level of OST1, and the inactivity of the 3 PP2C protein phosphatases (ABI1, ABI2, and HAB1).These elements form a positive feedback loop in which the enzyme activity of RBOH leads to sufficient ROS abundance, which inactivates all 3 PP2Cs, the inactivity of which is sufficient for activating the minimum functional level of OST1, which in turn activates RBOH.This positive feedback loop is self-sustaining, meaning that the corresponding active or inactive states can be maintained once established.We will refer to this stable motif as the Main Stable Motif  ), which is consistent with the CO 2 -independent high stomatal conductance of the ht1-3D and ht1-8D dominant mutants [12,13].Importantly, the MSM and the 3 conditionally stable motifs are mutually exclusive.Since the lock-in of a stable motif is irreversible, the choice between them constitutes the decision-making process of the CO 2 signaling network (Fig 4C) . The system has 1 additional stable motif and 3 additional conditionally stable motifs that contribute in specific circumstances; these are described in section II of S2 Text.Identifying the stable motifs allows comprehensive identification of all attractors of the system, using the algorithm we developed in [29].We first confirm the unique Closure ON attractor in the presence of high CO 2 (Fig 4D ), where high CO 2 leads to the locking-in of the MSM and precludes the activation of the 3 conditionally stable motifs.The MSM leads to the activation of anion flow through the SLAC1 channel and to K + efflux, as well as to the activation of Aquaporins, which together satisfy all requirements to activate the node Closure.In the case of ambient CO 2 , stable motif analysis identifies 3 attractors.Two of them have the node Closure OFF, corresponding to the absence of stomatal closure observed under ambient CO 2 .These 2 attractors differ only in the state of CPKs (either ON or OFF).Interestingly, the third attractor under ambient CO 2 has the node Closure ON and is highly similar to the closure ON attractor under high CO 2 -they differ only in the CO 2 in-component and SLAC1 CO 2 (see Fig 4D, row 2 versus row 1).Trajectories starting from the pre-stimulus initial state used in our simulation cannot reach this attractor under ambient CO 2 , consistent with biological reality and with the attractor found by simulation (Fig 3C).Yet, a different initial state, e.g., one in which HT1 is inactive, can lead to this attractor (see S1 Text).An intervention (e.g., sustained activation of a node) can also induce stomatal closure under ambient CO 2 , as we will present later.
Logic diagrams elucidate signal transduction mechanisms and the roles of driver nodes.To better understand how the system reaches its attractors (e.g., the closure ON attractor under high CO 2 ), we introduce logic diagrams illustrating the signaling mechanisms (Fig 5 ).A logic diagram (see Methods) is an integration of the network and node states that yield a specific outcome (e.g., the ON state of the node Closure), constructed based on the framework of the "expanded network" developed by the Albert group [86,93,94].Edges in the logic diagrams are logical implications (i.e., sufficient activations), and black dots represent AND gates.To make the logic diagrams easier to parse, we also applied network reduction techniques that do not change the dynamic behavior of the system [95], and we merged nodes that behave similarly or as a group (e.g., PP2Cs, see Methods).
As shown in Fig 5A, high CO 2 activates the MSM (shown with bold font), which leads to ion and water efflux, and thus to an attractor with the node Closure ON.Fig 5B shows the mechanism for the lack of stomatal closure, e.g., under ambient CO 2 .When HT1 is ON, the conditionally stable motifs (bold-font nodes, with their dependence on HT1 = ON shown with red dotted edges) can lock in and cause the inactivity of the SLAC1 and QUAC1 channels, thus reaching a Closure OFF attractor.The knockout of OST1 also activates the conditionally stable motifs and leads to a similar Closure OFF attractor (see the third entry of Fig 4D).In addition to attractors with constant closure values (constant ON or OFF), the system also allows attractors with an unstable closure value that oscillates between ON and OFF indefinitely.Fig 5C shows the logic diagram of the predicted mechanism of unstable Closure state that could arise from the sustained inactivity of the SLAC1 channels (due, e.g., to knockout of GHR1, see the fourth entry of Fig 4D).The oscillation of Ca 2+ cyt leads to oscillating QUAC1 and K + efflux, which in turn lead to an oscillating Closure state.Had SLAC1 been active, it would have sustained anion efflux, leading to sustained Closure ON.
In addition to finding attractors, stable motif analysis can also pinpoint key nodes called drivers, which determine the decisions in the system: once a driver node is sustained in a certain state, it can lock in a corresponding stable motif, forcing a decision in the system.Thus, external control of driver nodes can drive the system from one attractor to another [27,29].For example, ROS (marked in bold outline) is an internal driver of the MSM, as it inhibits all 3 PP2C proteins, which then allows OST1 activity, which in turn activates RbohD/F, closing the cycle of self-sustained activity.Furthermore, CaIM (Ca 2+ influx through the membrane) is an external driver of the MSM, as it inhibits the PP2Cs.When ROS or Ca 2+ is externally provided, it will lock in the MSM, thus lead to a closure ON attractor, even in ambient CO 2 (see Fig 4D).This explains the experimentally observed ROS or Ca 2+ -induced stomatal closure under ambient CO 2 [21,79,96].Moreover, ~HT1 (the OFF state of HT1) is also an external driver of the MSM, and leads to the same closure ON attractor as providing ROS or Ca 2+ (see Fig 4D), which explains the low stomatal conductance observed in the ht1-2 kinase dead mutant regardless of the CO 2 concentration [11].

Our model predicts how the perturbation of signaling elements affects high CO 2 response
A model is most helpful in predicting unknown behaviors of a system.We first demonstrate the predictive power of our model by simulating the guard cell response to high CO 2 under systematic single node perturbations.These simulations reflect how mutations or interventions could affect high CO 2 -induced closure.We simulate all single interventions that abolish a node's expression or activity as maintaining the node in the state 0 (OFF) and denote it "KO"; we also simulate all single interventions that make a node constitutively active by maintaining the node in the state 1 (ON) and denote it "CA."To quantify the stomatal closure response, we develop an improved approach from [3] combining 3 metrics: the final closure state, the cumulative percentage of closure (CPC), as well as the attractor repertoire of the system (see Methods).The CPC value represents the fraction of simulations in which Closure = 1, summed over 30 timesteps, thus it ranges from 0 to 30.As before [3,54], we observe and categorize 4 typical behaviors different than the wild type (WT) response: (1) we call the behavior in which the percentage of closure increases faster than in the WT system "hypersensitivity"; and (2) the behavior in which the percentage of closure increases more slowly than in the WT system "hyposensitivity"; (3) we refer to the behavior in which the percentage of closure converges to or oscillates around an intermediate value "reduced closure"; (4) and finally the category "no closure" is the lack of closure in simulated stomata.We present the summarized high CO 2 closure response categories under systematic perturbations in Fig 6A .The complete perturbation results can be found in S4 Table.We illustrate simulated response curves of these categories in Fig 6B : hypersensitivity under external Ca 2+ (CaIM CA), hyposensitivity under CaIM KO, reduced closure under H + -ATPase KO, and no closure under MPKs KO (mpk4/ mpk12 mutant).These perturbation response categories offer a way to characterize each signaling element's contribution to high CO 2 -induced closure, from weakest (nodes whose KO and CA are in the "Close to WT" and whose CA is in the "Hypersensitivity" category, such as ROS).
The 5 response categories are well-explained by our stable motif and network analysis.We find that in both hyper-and hypo-sensitivity cases, the final attractor of the system is almost the same as that of the wild-type system under high CO 2 ; the final attractor only differs in the state of the perturbed node and up to a few other nodes regulated by it.As illustrated in the logic diagrams (Fig 5) earlier, the typical reason for hyper-sensitivity (i.e., faster closure) is the immediate accomplishment of a node state that would naturally happen only after a previous event; e.g., CaIM normally would only happen after GHR1 activation (Fig 5A), but it happens immediately with CaIM CA.The mechanism for slower closure is the elimination of one of 2 redundant signaling pathways: for example, CaIM KO will eliminate one pathway to increase the cytosolic Ca 2+ level (a necessary condition of closure); the alternative pathway to activate Ca 2+ cyt is by CIS (Ca 2+ influx from stores).The model prediction that perturbation of specific anion channels, e.g., QUAC1 KO or R256 residue KO (interpreted as SLAC1 CO 2 = OFF) leads to hyposensitivity is due to the same type of mechanism.Reduced closure occurs when the state of the node Closure oscillates in all attractors, such as for GHR1 KO (see row 6 in Figs 4D and 5C).GHR1 KO leads to the deactivation of the SLAC1 anion channel, leaving QUAC1 as the only conduit of anion flow, and the activity of QUAC1 follows the Ca 2+ cyt oscillation, resulting in an attractor with oscillatory closure.The final category of behavior, "no closure," is the lack of closure in all simulated stomata; the system converges into an attractor (or mix of attractors) with closure OFF.The mechanism that leads to this outcome is disruption of ion or water efflux, either through a mutation in a transporter (e.g., Aquaporins KO in Our model predicts interventions that yield closure even under ambient CO 2 .We performed systematic perturbation simulations of the model under ambient CO 2 (i.e., with the node "high CO 2 " set to OFF) and found 18 interventions that can result in Closure = ON in the simulations (Fig 6C ; comprehensive results are in S4 Table ).Our stable motif analysis identifies that the reason for the effectiveness of these interventions is the activation of the MSM coupled with the inactivation of HT1 (see Fig 5A).HT1 KO itself is an example of such a stable-motif-activating intervention; this result agrees with the closed stomata or low stomatal conductance observed in the ht1-2 mutant [11,13].Thus, the model predicts that constitutive In agreement with the model predictions, it was observed experimentally that provision of ROS [79] or CaIM (Ca 2+ influx through the membrane, which is equivalent to providing external Ca 2+ ) [21] results in closure under ambient CO 2 .The predicted closure-inducing effect of additional interventions is supported by experimental evidence, e.g., NO [62,84], cADPR [85], 8-nitro-cGMP [38] were shown to induce closure as a signaling element of the ABA-induced closure network.

Experimental validation of model predictions
Our model predicts that the activation of the MSM is an important part of the process of stomatal closure both in ambient and high CO 2 (see Figs 5A and 6D).A follow-up conclusion from this prediction is that earlier activation of the MSM would lead to more rapid closure.One of the nodes whose constitutive activity is predicted to activate the MSM and lead to closure in ambient CO 2 is NO.In our model, applying NO under high CO 2 activates the MSM sooner.Hence, supplying NO is predicted to cause hypersensitivity of closure under high CO 2 (i.e., faster high CO 2 -induced closure).To assess the validity of this prediction, we experimentally tested the effect of NO on high CO 2 response by applying SNP (sodium nitroprusside dihydrate), an NO donor (Fig 7A, see Methodss for experimental assay description).To observe any acceleration of the closure process, we used short treatments (10 and 20 min).We found that under high (800 ppm) CO 2 , the stomatal aperture under SNP treatment (gray bars) is significantly lower than that without SNP treatment (blue bars, p = 2.8e-2 for 10 min and p = 3.4e-2 for 20 min, Student's t test).Furthermore, the mean aperture corresponding to the 10 min combined treatment of 0.15 mM SNP and high CO 2 is lower than the mean aperture corresponding to the 20 min high CO 2 treatment (rightmost blue bar, p = 4.0e-2, Student's t test), indicating that SNP treatment speeds up the high CO 2 -induced closure process, consistent with our model prediction.As further confirmation, we performed a similar experiment using the alternative NO donor SNAP (S-nitroso-N-acetylpenicillamine) and found similar results (see details in S2 Text section III).
Furthermore, we experimentally validate the prediction that NO can induce closure under ambient CO 2 .Our model predicts that under ambient CO 2 , NO can still activate the MSM, consequently inducing closure (Figs 4C and 6D).Previously, NO has been shown to induce closure as a signaling element of the ABA-induced closure network [62,84].Yet, it has been unknown whether the degree of NO-induced closure (aperture reduction) is comparable to that of high CO 2 -induced closure.We found that applying the NO donor SNP under ambient Taken together, these experiments not only validate our model predictions regarding the effects of NO, but also support the model-predicted importance of the MSM.

Stable motif analysis predicts a new regulatory relationship between NO and ABI2
Analysis of the inconsistencies between the model's results and experimental observations can identify edges that are missing from the network model.The node NO also offers an example of this kind.For example, while our model recapitulates NO-induced closure in ambient CO 2 , it does not capture the observation that NO depletion leads to loss of stomatal closure under elevated CO 2 in tomato [24].We performed this experiment in Arabidopsis and observed that 200 μM of the NO scavenger cPTIO impairs high CO 2 -induced stomatal closure (Fig 8A).Under high CO 2 , the aperture size with cPTIO treatment is significantly different than that without cPTIO treatment (p = 1.9e-3,Student's t test).This observation suggests that there is a not-yet-identified regulatory effect of NO that makes it necessary for a closure mechanism.This unidentified regulatory effect of NO can be predicted from our stable motif analysis.In our current model under high CO 2 and NO KO, the MSM can still lock-in by the activation of its driver, the inactivity of HT1, thus leading to a closure ON attractor (Figs 8B and 4D).According to current knowledge, NO participates in the pathway that leads to calcium influx from intracellular stores into the cytosol (CIS) [97].Disruption of this pathway eliminates one mechanism of increasing the cytosolic Ca 2+ level, but it does not perturb the MSM.To be able to drive the system away from a closure ON attractor, NO KO should prevent the lock-in of the MSM that activates closure (Fig 8C).Thus, we hypothesize that NO regulates at least one of the nodes in the MSM, namely, it activates ROS, RbohD/F, or OST1 minimum function, or inhibits ABI1, ABI2, or HAB1.Indeed, we experimentally found that NO inhibits ABI2 phosphatase activity (Fig 8D and see Methods).This result suggests that the mechanism (or one of multiple contributing mechanisms) leading to impairment of stomatal closure under NO knockout could be the increase of ABI2 activity, which obstructs the minimal functional OST1 activity and precludes the lock-in of the MSM.This newly discovered inhibitory regulation likely contributes to NO-induced stomatal closure by the NO-induced decrease in the phosphatase activity of ABI2 making the activation of the MSM easier.The experimental confirmation of this stable-motif-based model prediction highlights the power of our model in predicting/ prioritizing potential regulatory relationships and revealing the mechanisms of signal transduction.

Our model predicts which nodes of the high CO 2 signaling network could be potential targets of currently unconnected nodes
There are signaling elements that are convincingly implicated in CO 2 signaling, yet their connection to the signaling network remains unknown.During our literature review and data curation process, we placed such elements into the "insufficient connection" category.For example, big1 mutants are compromised in elevated CO 2 -induced stomatal closure and bicarbonate activation of S-type anion channel currents [98], yet BIG1's connection to known guard cell signaling elements remains undiscovered.Connecting these elements to the existing CO 2 network can improve understanding of guard cell CO 2 signaling and could lead to testable predictions.Our model offers a way to identify and prioritize potential targets of these unconnected elements based on the effects of node perturbations (Fig 6).The idea is that, if a node perturbation yields a similar response in the model to the observed effect of the unconnected signaling element, the node is then a potential target of the signaling element.For example, the big1 mutant is known to have impaired CO 2 -induced closure.Our model simulations (Fig 6A) yield a comprehensive list of 24 single node perturbations that could impair closure.We predict that the big1 mutant regulates one of these nodes; any alternative single target of big1 would not cause loss of closure in the model.Moreover, one can utilize additional evidence to further narrow down the target candidates.For example, as big1 mutants are also compromised in bicarbonate activation of S-type anion channel currents, we can further filter the 24 candidate nodes to those that satisfy both conditions: (1) the node is downstream of bicarbonate; and (2) the node's perturbation results in reduced SLAC1 activity (reflected by the node AnionEM, which incorporates the effect of both SLAC1 channel nodes).This filter narrows the list of potential targets of BIG1 to 9 nodes, namely, MPKs, ABI1, ABI2, GHR1, HT1, OST1 minimum function, RbohD/F, and ROS.One can design experiments to test if some of these nodes are indeed perturbed in big1 mutants.If a broader pool of candidate targets is needed, one could also consider the nodes whose perturbations causes hyposensitivity in high CO 2 -induced closure (5 nodes in the last row of Table 4).These nodes may be appropriate because it might be difficult to distinguish reduced sensitivity and hyposensitivity in experiments.
By a similar logic, the "hypersensitivity" response category in Fig 6A allows target prediction for signaling elements whose perturbation (knockout or constitutive activation) promotes closure.For example, faster closure is observed in plants lacking the ABC transporter AtABCB14 [99].Hence, our model predicts that the knockout of AtABCB14 causes a node perturbation that leads to hypersensitivity, as annotated in Fig 6A .Furthermore, one can predict/prioritize targets of a signaling element as those that can induce closure in ambient CO 2 , based on the node perturbations that can include closure in ambient CO 2 (Fig 6C).To summarize, our model provides ways to prioritize potential targets to link previously unconnected nodes with the high CO 2 signaling network.Importantly, we only used BIG1 and AtABCB14 as examples; the proposed prediction method is general.

Discussion
In this work, we constructed a signaling network that integrates current information regarding the signaling mechanisms underlying high CO 2 -induced stomatal closure in guard cells.We assembled the nodes and edges of the network model from a thorough literature review process, with rigorous definitions of which signaling elements are involved in CO 2 signaling and which are not; we also deployed a multi-node presentation to capture multiple mechanisms of the same element.Finally, we conducted a comprehensive model performance evaluation, determining that our model indeed has high consistency with experimental observations.A previous model by Schroeder and colleagues of ABA and high CO 2 -induced closure, constructed as a case study of the Boolink graphical user interface [31], added to a model we originally created for ABA-induced stomatal closure [3] 6 nodes (CO 2 , CA1, CA4, MPK4/12, HT1, and CBC1/CBC2) and changed 5 functions (GHR1, CaIM, Ca 2+ c , Microtubule, and H 2 O efflux) without comprehensive evaluation of the model's accuracy in recapitulating all experimental observations regarding high CO 2 -induced closure.In fact, 2 of the additions make the rest of the network superfluous: specifically, the function that controls CaIM has "OR CO 2 " in it, which means that high CO 2 induces CaIM regardless of the rest of the regulators of CaIM; the function of H 2 O efflux also has "OR CO 2 " in it, which makes water efflux happen directly in response to CO 2 regardless of the rest of the regulators of aquaporins.The combined effect of these 2 assumptions results in closure.In other words, the model of Karanam and colleagues short-circuits essentially the entire CO 2 signaling network by making CO 2 directly activate the very downstream part of the network.As a result, the model will not be able to capture the experimentally observed responses to many perturbations in the signaling network, e.g., the impaired high CO 2induced closure under slac1 [81], ht1-2 [12], and mpk4/mpk12 [14] mutants.By comparison, a rigorous model construction process allows our CO 2 model to comprehensively reproduce experimental observations related to high CO 2 -induced closure.Table 4. Model-predicted potential regulation targets to impair high CO 2 -induced closure.If a knockout of a signaling element impairs high CO 2 -induced closure, the model predicts the knockout will perturb one of these nodes.The first 3 rows are categorized based on the network location of the node, i.e., which component it is in.The last row is an extension where we hypothesize that perturbations causing hyposensitivity in the model could also be considered as having the potential to impair closure.Bold font indicates an example of 9 prioritized targets for BIG1, which are downstream of bicarbonate and whose perturbation can also cause reduced SLAC1 activity, as observed in the big1 mutant.

Network methodologies reveal key decision-making motifs, pinpoint driver nodes, elucidate signal transduction mechanisms, predict new edges and responses
Our network-based modeling framework integrates piece-wise (node to node) biological knowledge to form a comprehensive view of the signal transduction network underlying high CO 2 -induced stomatal closure.This framework complements conventional biological approaches that often focus on single regulators or single pathways.Our analysis supports recent findings that HT1 is a key decision-making node of the CO 2 signaling system [15,16] and offers additional insights regarding the ambient CO 2 setting: closure induced by, e.g., exogenous ROS or Ca 2+ in ambient CO 2 must involve the inactivation of HT1 and the activation of a positive feedback loop (the MSM, formed by ROS, OST1, and the PP2Cs).Conversely, under conditions that activate HT1 the nodes of this positive feedback loop may also lock-in to their opposite states and lead to closure OFF.
We have also introduced the concept of logic diagrams to visually depict the logical implications of events in the closure process, highlighting key nodes and pathways (subgraphs) responsible for decision-making of the system.The logic diagrams explain not only how signal transduction leads to stomatal closure under high CO 2 , but also how closure is lost or restored.Leveraging these concepts, we explain the mechanisms behind the 5 types of simulated closure responses in the presence of perturbations.
Combining these network-based methodologies allows the generation of model predictions that lead to novel experimental findings.It is important to note that inconsistencies between the model and experimental results also yield insights.For example, we found that the model is incapable of reproducing the experimental observation that NO knockout impairs high CO 2 -induced closure; by understanding that this inconsistency results from the NO-independent lock-in of the MSM, we predicted that NO must regulate the MSM, and successfully validated this prediction experimentally, showing that NO inhibits ABI2, one of the PP2Cs in the MSM.

Unexplained mechanisms in CO 2 signaling
There are signaling elements that are known to mediate closure in response to a certain signal, but are understudied in CO 2 signaling.In our model construction process, when we designate each piece of knowledge regarding the regulation of a node by another into one of the 5 categories (see Tables 2 and S1), we place such elements into the "generic" or "insufficient evidence" categories.This way, our categorization of evidence can prioritize experiments to explore the roles of signaling elements in the context of CO 2 signaling.
Past studies on the convergence of the signaling networks corresponding to ABA-induced stomatal closure and high CO 2 -induced stomatal closure yielded inconsistent results.This convergence was first suggested by Webb and Hetherington [5].Xue and colleagues [9] found that the quadruple ABA receptor mutant pyr1;pyl1;pyl2;pyl4 showed slower closure but no other impairment under 800 ppm CO 2 .On the contrary, Chater and colleagues [22] and Dittrich and colleagues [100] found that quintuple mutants of ABA receptors from the PYR/PYL/ RCAR family exhibit impaired high CO 2 -induced closure, and that ABA biosynthesis-deficient nced3;nced5 plants were impaired in high CO 2 -induced closure.Movahedi and colleagues [101] reported hypersensitivity to CO 2 in cyp707a1 cyp707a3 ABA-hyper-accumulating mutants, supporting the convergence of the 2 pathways.By contrast, there is also evidence for ABA-independence of high CO 2 signaling [2], e.g., OST1/SnRK2.6 protein kinase activity in guard cells is not activated by CO 2 elevation [34,35], and cytosolic pH does not increase under high CO 2 [9,32].Moreover, CO 2 regulation of PYR/PYL/RCAR signaling activity, as well as a mechanism by which PYR/PYL/RCAR receptors would affect high CO 2 -induced closure independently of changes in ABA concentration is not known.Considering all the above evidence, we placed RCARs in the "insufficient connection to CO 2 " category when constructing this model.Nevertheless, our model can be expanded to include a node that represents endogenous ABA if one assumes that it is necessary for OST1 minimum function but not sufficient to inhibit the PP2Cs and we now discuss this scenario in S2 Text section IV.
Our model construction and validation process included the examination of multiple different hypotheses regarding the role of certain nodes in high CO 2 signaling.For example, the role of RHC1 in CO 2 signaling is under debate: while Tian and colleagues showed experimental evidence that RHC1 is an essential CO 2 signaling component [17], Tõldsepp and colleagues observed that a different rhc1 mutant shows no reduction of high CO 2 -induced closure.We tested both versions of the HT1 regulatory function, with and without RHC1, and found that the function without RHC1 has a higher accuracy in capturing experimental observations.Therefore, we adopt the HT1 function without RHC1, more consistent with [18].Our interpretation is that RHC1's regulatory effect on HT1 is weaker than that of MPKs, thus it is insufficient to include in a Boolean model.Some CO 2 signaling mechanisms remain unclear.For example, observations about guard cell cytosolic pH suggests undiscovered buffering capacity in the cytosol.Xue and colleagues and Brearley and colleagues [9,32] observed that cytosolic pH does not increase under high CO 2 -induced closure.On the other hand, after CO 2 enters the guard cell, it undergoes the chemical reaction CO 2 + H 2 O !HCO 3 -+ H + catalyzed by CAs, producing H + .So, there must be some buffering capacity or other mechanism to keep cytosolic pH constant.Understanding this mechanism may affect other components in the network, and its refinement may improve the model as well as our understanding of the high CO 2 signaling mechanism.Future discovery of new signaling elements and interactions that participate in high CO 2induced closure will call for updates of this model, including the addition of new nodes and/or new edges, the possible removal of assumed edges (shown as dotted lines in Fig 2), or changes to the regulatory functions of certain nodes.It is important to note that even when new signaling elements and interactions are incorporated into the CO 2 model, certain structural and dynamic features of the model will be preserved.For example, the network will still have a strongly connected component (see Fig 1B) that includes at least the 20 nodes currently in it.The updates will be such that the model preserves the current model's agreements with experimental results.The coupling of the structural and dynamical preservation means that the decision-making role of the MSM will also be maintained in the updated model, even if newly added edges will increase its size (for example, by NO joining it).It is likely that newly discovered interactions will support the regulatory relationships that drive the decisions of the system.For example, the inhibition of ABI2 by NO is an additional mechanism by which NO contributes to the activation of the MSM.Ultimately, a cycle between experimental and modeling approaches offers the most effective way to elucidate the full network of high CO 2 -induced stomatal closure.

Construction of the CO 2 -signaling network and dynamic model
As described in the main text, we extract regulatory relationships from the literature (S1 Table ) and compile them into nodes and edges of the network (S2 Table ).Then, we construct the regulatory functions (S1 Text) based on the regulatory relationships.

Simulation settings of the dynamical model
In simulation, we initialize the in silico stomata in a pre-stimulus open state, with closureinhibiting nodes ON (namely ABI1, ABI2, HAB1, HT1, and H + ATPase) and other nodes OFF (see S1 Text for complete initial condition).For example, Anion Efflux, a node that becomes activated during high CO 2 -induced closure, is initialized in the OFF state.We deploy a discrete-time simulation, with stochastic random order asynchronous update: at each discrete round of updates (or timestep), all nodes are updated in a randomly selected order.A different order is randomly generated in each timestep.The main advantage of random order asynchronous update is that it can capture consensus dynamics without implementation of exact timing, as the actual timing and durations of the processes happening in the system are typically unknown.The stochasticity of the update helps avoid spurious behaviors that depend on synchrony (i.e., which only appear under synchronous update, where all nodes are updated at the same time, which is biologically unrealistic).We perform this simulation for 30 timesteps, long enough for the system to reach its attractor (stable long-time behavior).In this way, we can simulate how a node's state changes over time during the stomatal closure process.
We perform 1,000 replicate simulations, i.e., 1,000 "in silico stomata" simulations, then calculate the "percentage of node activation" at each timestep as our metric, meaning the percentage of simulations in which the node (e.g., Closure) is in state 1 (ON).For example, if 1,000 replicate simulations yield 510 cases where closure is ON at timestep 30, the simulated percentage of closure is 51%.In this way, oscillations and multistability with different node states will display a percentage of node activation value between 0% and 100%.A sample of our simulation code is provided on GitHub: https://github.com/jack802gan/CO2-model. We also recommend Cubewalker [102] as a more recent and more efficient simulation tool for Boolean dynamic models.

Stable motif and attractor analysis
Stable motifs are the smallest strongly connected components in a system that can sustain a unique steady state for the constituent nodes [28,86].Stable motifs are identified in an expanded representation of the network (which includes the network's regulatory functions as part of the network structure).A weaker version of a stable motif, called a conditionally stable motif [92], can sustain a unique steady state for its nodes as long as certain node(s) that regulate the motif are in a specific fixed state.The stable motifs and conditionally stable motifs of a Boolean network determine its attractors: one can uniquely associate sequences of stable motifs (stabilized in the order given by the sequence) to each attractor.
A driver (or driver set) of a stable motif is the node (or node combination) and its state that, when kept fixed, eventually leads to the locking in of the stable motif after sufficient updates.We used the python library pystablemotifs [29] to analyze the stable motifs, attractors, and their control methods.The library is available through the GitHub page https://github.com/jcrozum/pystablemotifs.

Validation against experimental observations
We use the experimental observations in S1 Table that were not used in model construction (a subset of the indirect regulations, marked as "V" for validation in S1 Table ) as the ground truth for model validation.Because quantifying the (reduction of) stomatal aperture is difficult in both experiment and Boolean model, when comparing the experiment-based expected outcomes with the model simulation results, we employ a lenient criterion for consistency: any model-predicted defect (i.e., node percentage < 100%) is considered qualitatively consistent with either reduced response or no response.For example, if a simulation yields a closure percentage = 0.5, then it is considered consistent with both "reduced closure" and "no closure" in experimental observation.This type of consistency criterion was used in previous modeling as well [3].We provide the comprehensive validation results in S3 Table .The supplementary file also includes the trivially consistent cases due to circular reasoning, denoted "true by construction."For example, bicarbonate-induced activation of S-type anion channels was reduced in the dominant active abi1 mutants (ABI1 CA) [78], but because ABI1 is an input in both SLAC1 functions and ABI1 ON is sufficient to deactivate both SLAC1 nodes, this experimental observation is captured by model/rule construction, and calling it a validation would be circular reasoning.This type of validation is not counted in our assessment of model accuracy performance.

Definition of the response categories in CPC values, for simulated knockout or constitutive activity
We defined the CPC as the sum over 30 timesteps of the fraction of simulations that have the closure node in state 1 (ON), over 1,000 replicate simulations.Thus, CPC ranges from 0 (if the percentage of closure were 0% at every timestep) to 30 (if the percentage of closure were 100% at every timestep).The CPC metric helps capture time-sensitive behavior such as hypo-/ hyper-sensitivity, reflected as a small CPC difference from wild-type CPC.
We define the 5 categories of responses based on the CPC value and attractor state: close to wild type, hypersensitivity, hyposensitivity, reduced closure, and no closure.We define the wild-type response range from the default high CO 2 -induced closure response.Since the simulation is stochastic (i.e., the CPC value can be different for each simulation), we determine the default CPC range by 600 batches of 1,000 replicate simulations under the default high-CO 2 signal, and define its CPC range as the "close to wild type" CPC range (20.405 to 20.709).CPC values above this range are classified as hypersensitivity, meaning that the stomata reach the closure state significantly faster than wild type.CPC values below the wild type range with 100% percentage closure at the final timestep are classified as hyposensitivity.Reduced closure is the CPC value being much lower than wild type range, with the final percentage of closure being between 0 and 1.In all of the observed cases in this last category, the closure value is oscillating in the attractor.Lastly, no closure is defined as the CPC value being zero or very close to zero, with final percentage of closure being zero.A small positive CPC can occur here, because closure may temporarily turn ON during the response process, though it cannot be maintained.The comprehensive CPC simulation results are provided in S4 Table.

Logic diagrams of the CO 2 network model
The logic diagram is an integration of the network and node states that yields a specific outcome (e.g., the ON state of the node Closure).It is an expansion of the concepts of elementary signaling mode [93] and logical domain of influence [94], both of which are subgraphs of the expanded network (also called logical hypergraph), which incorporates the Boolean functions of all the nodes in the network [86].The logic diagram indicates the node states visually: nodes whose name is preceded by ~and have gray background are inactive, nodes with green background are active, and nodes whose name are followed by "osc" and have a white background oscillate.A single edge in a logic diagram represents a sufficient regulatory relationship that can achieve the indicated state of the target node; multiple edges incident on the same black dot represents "AND" logic, indicating that all nodes are required to achieve the target node state.A logic diagram may condense fully linear pathways of nodes and edges in the original network, to simplify visualization.This type of condensation does not change the regulatory logic of the network model [95].
To help visual parsing in the CO 2 network model, we also merged a few nodes of the original network.The node "PP2Cs" merges 3 members of the PP2C protein family, because these nodes have identical regulation.The technical definition of the merged node is PP2Cs = ABI1 or ABI2 or HAB1.The node "Ion flow" merges AnionEM and K + efflux (i.e., Ion flow = AnionEM and K+ efflux).We also use a single node "Ca 2+ cyt osc" to represent the Ca 2 + cyt -Ca 2+ ATPase negative feedback cycle, capturing the oscillating behavior of Ca 2+ cyt in the presence of a Ca 2+ source (CaIM or CIS) and due to the inhibitory effect of the Ca 2+ ATPase.We use hexagonal symbols for "PP2Cs" and "Ion efflux" as they represent a node set (Fig 5).Note the "Ion efflux" hexagon node needs multiple incident edges to activate, causing an exception of the edge definition.
The logic diagram can be read as a sequence of events and the conditions necessary for these events.For example, the sequence of events necessary for high CO 2 -induced closure described by Fig 5A is the following: high CO 2 activates HCO 3 -, which inactivates HT1 via the formation of the MPKs-HT1 complex.HT1 being OFF then activates OST1 minimum function, which eventually activates the MSM.The simultaneous activation of the MSM and inactivation of HT1 (indicated by the black dot) activates GHR1, which in turn activates CPKs.GHR1 activates CaIM, which leads to cytosolic Ca 2+ transients or oscillations.The MSM, GHR1, Ca 2+ cyt , HCO 3 -, as well as the inactivity of CBC1/2 contribute to the activation of Ion flows.The simultaneous activation of Ion flows and OST1-activated Aquaporin yield H 2 O efflux, which then leads to Closure.As another example, the logic diagram in Fig 6D shows that the activation of the MSM (e.g., by CaIM) leads to CIS, which then leads to cytosolic Ca 2+ transients or oscillations, which activate MPK activity, which inhibits HT1, which then allows the activation of GHR1, then of the generic SLAC1 node, and finally the activation of the node Closure.

Plant materials and growth conditions
Arabidopsis thaliana Columbia (Col-0) accession plants were used as wild type.Col-0 were germinated on MS agar plates for 2 weeks and then grown in soil (1:1 mixture of Metro Mix 360: Sunshine Mix LC1 (Sun Gro Horticulture) in a growth chamber with an 8-h-light/16-hdark cycle, 150 μmol photons m −2 s −1 light, and 22˚C (day)/19˚C (night)).Young fully expanded leaves from 4-to 5-week-old plants were used for stomatal aperture assay experiments.

Stomatal aperture assay
Stomatal apertures were assayed as described previously with slight modifications [103].In brief, epidermal peels were hand-peeled from the abaxial surfaces of leaves and then were floated interior-side down on an opening solution (20 mM KCl, 10 mM MES-KOH (pH 6.15), 50 μM CaCl 2 ) under white light (150 μmol m −2 sec −1 ) for 3 h to induce maximum opening of the stomata (initial aperture).After the 3 h preincubation period, epidermal peels were floated in a six-well plate (each well with a diameter of 35 mm and height of 17 mm) containing preequilibrated high (800 ppm) or ambient (400 ppm) CO 2 opening solution (approximately 3 ml).Subsequently, the multiwell cell culture plate was placed in a petri dish chamber and sealed with parafilm.The high or ambient CO 2 concentration was continuously maintained within the chamber using a Licor-6400 (flow rate of 500 μmol s −1 ) while the samples were further incubated under white light (150 μmol m −2 sec −1 ).After 3 h of incubation, the epidermal peels were imaged with a 40× objective using a light microscope (Nikon Diaphot 300) coupled to a digital camera (Nikon Coolpix 990).Stomatal apertures were measured in each image using ImageJ.To determine the effects of NO donor SNP (sodium nitroprusside dihydrate; Sigma, catalog number 71778), SNAP (S-Nitroso-N-acetyl-DL-penicillamine; Sigma, catalog number N3398), or NO scavenger (2-(4-carboxyphenyl)-4,4,5,5-tetramethylimidazoline-1-oxyl-3-oxide; Sigma, catalog number C221) on CO 2 -induced stomatal closure, experiments were performed according to the method described by [104] with some modifications.Briefly, epidermal peels were incubated for 3 h in the opening buffer to induce stomatal opening.Subsequently, the epidermal peels were transferred to pre-equilibrated high (800 ppm) or ambient (400 ppm) CO 2 opening solution (approximately 3 ml) added without or with NO donor or NO scavenger at the indicated concentrations.Then, the epidermal peels were treated under ambient CO 2 (400 ppm CO 2 ) or elevated CO 2 (800 ppm CO 2 ) in light for the indicated durations.A total of at least 30 stomatal apertures were measured per treatment for each individual experiment.The datapoint for each individual experiment is the mean aperture of these stomata.Data presented in Figs 7A, 7B, and 8A are the mean and standard error in the mean of the 3 data points coming from 3 independent experiments.The relevant data are provided in S1 Data.

Phosphopeptide-based phosphatase activity assay
To examine the effect of NO donor on ABI2 phosphatase activity, GST-tagged ABI2 was heterologously expressed in Escherichia coli and subsequently purified as described [105] with minor modifications.Briefly, pDEST15-ABI2 was constructed by transferring ABI2 cDNA from pCR8/GW/TOPO TA vector to the N-terminal GST fusion expression vector pDEST15 (Invitrogen) by Gateway LR recombination (Invitrogen).The cloned construct was verified by sequencing.The GST-ABI2 protein was expressed in E. coli BL21(DE3) pLysS Rosetta cells by induction with 0.5 mM IPTG for 16 h at 16˚C.The cells were harvested by centrifugation (6,000 g, 8 min, 4˚C) and stored at −80˚C until use.All purification procedures were carried out at 4˚C.The final cell pellet was resuspended in 10 ml of B-PER reagent (Pierce, catalog number 78243) supplemented with 1 mg/ml lysozyme (Sigma, catalog number L6876), 25 mg/ ml DNAse and EDTA-free complete protease inhibitor (Thermo Fisher Scientific, catalog number A32955), and incubated for 60 min at 4˚C with gentle shaking.The cell lysate was spun at 16,000g for 15 min at 4˚C and the resulting supernatant was filtered using a 0.45 μm filter.Subsequently, the supernatant was loaded onto a 5 ml QIAGEN polypropylene column packed with GST glutathione agarose resin (Pierce, catalog number 16100) and the columnbound protein was eluted with equilibration buffer (25 mM Tris-HCl (pH 7.5), 150 mM NaCl) supplemented with 20 mM reduced L-glutathione (Sigma, catalog number G4251).The purity of the recombinant GST-ABI2 was estimated by analysis of aliquots separated on a 12% SDS polyacrylamide gel stained with Gel-Code Blue (Thermo), and the concentration was estimated based on known concentrations of Fraction V BSA (Pierce, catalog number 23209) run on the same gel.
The effect of NO donor on ABI2 phosphatase activity was assayed using the Ser/Thr phosphatase assay kit (Promega; catalog number V2460).The reaction was performed in a 50 μl volume containing 28 nM of purified recombinant GST-ABI2 protein in the presence or absence of the indicated concentration of NO donor SNP.After 30 min incubation at 30˚C, the reaction was stopped by addition of 50 μl of molybdate dye/additive solution.The resulting mix was incubated for another 20 min at room temperature followed by absorbance measurements at 630 nm in a Synergy Neo2 multimode plate reader (Biotek).The amount of released phosphate was calculated based on a standard curve obtained using phosphate standard solutions of known concentration.The relevant data are provided in S1 Data.

Statistical testing
We used two-tailed Student's t tests when comparing the stomatal aperture measurements, to evaluate if a significant change of aperture size or phosphatase activity is observed in an experiment.The statistical tests were performed in Python (Anaconda JupyterLab).Relevant data and an example of code are provided on GitHub: https://github.com/jack802gan/CO2-model. DOI: https://zenodo.org/doi/10.5281/zenodo.10665442

Fig 1 .
Fig 1. Flow chart of the process of network construction and analysis.We construct a signaling network based on literature review, resulting in a network of 47 nodes and 95 edges.We then construct a Boolean dynamic model, which adds node states and regulatory logic to the network.We assess the accuracy of our network model by comparing its results with known experimental observations, finding that the model simulations successfully recapitulate 91% of experimental observations.We also perform multiple network-based analyses on the model.Finally, we experimentally validate selected model predictions.

Fig 2 .
Fig 2. (A) The high CO 2 signaling network with 47 nodes and 95 edges.Yellow nodes represent validated elements in CO 2 signaling, i.e., category (1); white nodes represent assumed elements based on generic evidence, i.e., category (2).(B) Acyclic graph representation of the network highlighting the in-component, the strongly connected-component and the out-component, forming a hierarchical structure of signal transduction.(C) Illustration of the multi-node representation of signaling elements that participate in multiple pathways and mechanisms, using MPKs as an example.MPKs are known to inhibit HT1 by a physical interaction induced by bicarbonate and independent of the MPKs' kinase activity[15,16], and also are necessary in Ca 2+ -induced SLAC1 activation[37].
Fig 3B), consistent with the expectation that stomata are open in the absence of CO 2 elevation.The row marked "ambient CO 2 " in Fig 3C shows the states of all nodes in the attractor under ambient CO 2 .

Fig 3 .
Fig 3. (A) Example demonstrating how a Boolean regulatory function is constructed from observations.(B) A demonstration of model simulation, capturing stomatal closure under high CO 2 , and lack of closure under ambient CO 2 .(C) Summary of node states in the attractor with closure ON under high CO 2 , and the attractor with closure OFF under ambient CO 2 , found from the model simulations.A green box represents the ON state, gray represents the OFF state, and X represents oscillation.If a prior experiment compared the states of a node under high and ambient CO 2 , we mark the corresponding cells with "e."In all such cases, the model recapitulates the experimental result.https://doi.org/10.1371/journal.pbio.3002592.g003

Fig 4 .
Fig 4. Stable motifs of the network model.(A) The MSM that leads to closure in the dynamic model.(B) Superposition of the 3 overlapping conditionally stable motifs that lead to attractors that represent open in silico stomata.Each conditionally stable motif contains a different PP2C protein family member, indicated by dashed, dashdotted or dotted outline and edges, respectively.In the stable motifs, nodes whose name is preceded by ~are inactive in the corresponding motif.The black circles indicate that the regulatory effects incident on a node form an AND gate; for example, ABI1 will activate only if ROS and CaIM are simultaneously inactive.The red edges starting from HT1 indicate that HT1 = ON is the condition necessary for the conditionally stable motifs.(C) Decision-making diagram of the CO 2 signaling network.(D) Summary of the node states in the attractors of the system.Entries that contain a white box with * summarize multiple attractors in which the marked nodes take either on or off states.MSM, Main Stable Motif; ROS, reactive oxygen species.https://doi.org/10.1371/journal.pbio.3002592.g004

Fig 5 .
Fig 5. Logic diagrams show the mechanisms that lead to various stomatal behaviors.(A) Logic diagram of high CO 2 -induced stomatal closure.It also indicates that a closure attractor can be reached under ambient CO 2 if the MSM is activated by an internal or external driver.Three such drivers (out of 17) are indicated by bold outlines.(B) Diagram of the mechanism for the lack of closure due to the activation of the conditionally stable motifs, e.g., under ambient CO 2 or for ost1 knockout.(C) Logic diagram of the predicted mechanism of unstable closure that would arise with inactivity of the SLAC1 channels.Green nodes indicate the active state and gray means the inactive state; white nodes oscillate.MSM, Main Stable Motif.https://doi.org/10.1371/journal.pbio.3002592.g005 Fig 4D), or indirectly through activation of the non-closure conditionally stable motifs (e.g., OST1 KO in Fig 4D; the mechanism is illustrated in Fig 5B).

Fig 6 .
Fig 6.Systematic perturbation under high CO 2 with illustrations.(A) Response categorization for systematic perturbations of all nodes, both knockout (KO) and constitutive activation (CA) perturbations.(B) Simulation curves under high CO 2 that exemplify response types: normal WT closure, hypersensitivity under external Ca 2+ (Ca 2+ influx across the membrane (CaIM) CA), hyposensitivity under CaIM KO, reduced CO 2 sensitivity for H + -ATPase, or insensitivity to CO 2 under MPKs KO. (C) Model simulation predicts 18 interventions that could induce stomatal closure under ambient CO 2 .Some of these interventions have been validated to induce closure under ambient CO 2 , and some interventions have supporting evidence in previous publications.(D) Logic diagram of the mechanism of closure in ambient CO 2 .The node "MSM" is a merged representation of the MSM.Bold outlines indicate 2 of the 18 interventions that can simultaneously activate the MSM and inactivate HT1 under ambient CO 2 .For example, NO leads to oscillations in Ca 2+ cyt , which lead to MPK activity, which inactivates HT1, which then activates the MSM and leads to closure.MSM, Main Stable Motif; NO, nitric oxide.https://doi.org/10.1371/journal.pbio.3002592.g006 CO 2 can indeed induce stomatal closure (Fig 7B, gray bars), as predicted by our model.The reduction in aperture caused by 200 μM SNP is statistically significant (p = 1.6e-3,Student's t test) and close to the reduction in aperture caused by CO 2 of 800 ppm (black bars).This result confirms the model prediction that NO can induce closure in ambient CO 2 .

Fig 7 .Fig 8 .
Fig 7. Experimental validation of the predicted closure responses when applying NO (via the NO donor sodium nitroprusside dihydrate, SNP) in ambient or high CO 2 .(A) As predicted by the model, applying SNP causes hypersensitivity in high CO 2 -induced stomatal closure: the combined treatment (orange and gray bars in the 800 ppm CO 2 condition) leads to a lower mean aperture than high CO 2 alone (blue bars), and a 10 min combined treatment leads to a lower mean aperture than a 20 min treatment of high CO 2 alone.(B) Applying SNP (NO donor) under ambient (400 ppm) CO 2 induces closure in a dose-dependent manner (gray bars; 400 ppm CO 2 ), as predicted.The highest concentration of SNP yields a similar degree of closure as high (800 ppm) CO 2 (black bars).The bars indicate the mean ± SE aperture from at least 3 independent experiments with at least 30 apertures measured per treatment.Significance of aperture differences is determined by Student's t test: ** P < 0.01, * P < 0.05.The data underlying the graphs shown in the figure are provided in S1 Data.https://doi.org/10.1371/journal.pbio.3002592.g007

Table 1 . Edge list of the high CO 2 -induced closure network. The
edges represent interactions and regulatory effects.The start node (first column) indicates the source of the interaction and the end node (second column) indicates the target of the interaction."Type" records the sign and type of the edge.Positive edge sign, "+", indicates positive/activating regulation and "-" indicates negative/inhibitory regulation.For edge type, we use "D" for direct regulation, "I" for indirect regulation, and "A" for assumed regulation.In rare cases, we use "D/A" to represent an edge that is supported by a direct interaction but has assumed properties, e.g., assumed direction or sign."Ref."column records the reference that supports the edge.Edges representing known biophysical/biochemical processes, e.g., AnionEM !Depolarization, do not have a reference.The table follows alphabetic order of the start node.

Table 2 . Signaling elements that differ in the networks that lead to closure in response to ABA versus high CO 2 . Category Nodes/signaling elements
[3]e full names of the ABA signaling elements can be found in[3].** Details and references for these elements are provided in S1 Table https://doi.org/10.1371/journal.pbio.3002592.t002

Table 3 . Comprehensive validation comparing simulations against known experimental observations.
In the 47 total simulations in the table, 43 are consistent and, 4 are inconsistent (marked with bold, underlined font), with accuracy 43/47 = 91%.

Table 3 .
[92]tinued)The lock-in of the MSM is observed in the closure ON attractor (the first row of Fig 3C).The opposite states of the nodes involved in the MSM, namely the inactivity of RbohD/F, ROS, the lack of OST1 function, and the activation of one of the PP2C protein phosphatases, when coupled with the inactivity of CaIM and GHR1, form 3 instances of a weaker type of stable motif, called a conditionally stable motif[92](see Methods).Fig4Bindicates the superposition of the 3 conditionally stable motifs and shows that they can lock-in if HT1 is ON, which in the wild-type system only happens under ambient CO 2 .The lock-in of these conditionally stable motifs leads to the attractor with closure OFF (second row in Fig 3C